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Abstract — Zusammenfassung 



Algorithm 48. A Fast Algorithm for Clusterwise Linear Regression. A fast implementation of a formerly 
[5] published algorithm. is given. 

AMS Subject Classifications: 62J05, 62H30, 65F20, 65D10-. 
Key words: Cluster analysis, linear regression. 

Algorithmus 48. Ein schneUer Algorithmus zur klassenweise linearen Regression. Fur einen fruher 
publizierten Algorithmus [5] wird eine schnelle Implementation angegeben. 



Let be given m observations {y iy a ik ) (i=l,...,vn, k=l, ...,/) with m>l Then this 
algorithm tries to find a partition C x , . . ., C n of given length n for these observations, 
i.e. C i cM={l,,.,4 \Cj\>l, CjnC k = $ for j^k, C r u... uC n = M, and re- 
gression coefficients x 01 ^*^, r->>w) (/ = 1, ■ such that 



This objective is reasonable when the number m of observations is relatively large as 
against to the number I of variables and/or when the observations might stem from 
different groups. For /=1 and a ix ^\ (z = l, ... m) you will have the well-known 
minimum variance criterion from cluster analysis [6] in one dimension. 

The above idea was considered in [5] and an inefficient program was given, too, that 
additionally had a small mistake [7] but could easily be extended for L p norms. The 
purpose of the present paper is to give a really efficient implementation of precisely 
the same algorithm. 



For larger values of m an exact optimum of the objective function cannot be found 
within reasonable computing times [6], The following heuristic method is nearly 
identical to that exchange method that is successfully used for the minimum 
variance criterion and the quadratic assignment problem. 



1. Problem and Purpose 




2. Numerical Method 
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Step 1: Choose some initial partition C l5 . . C n that is feasible, i.e. | C, | > /, and some 
starting observation i = i a . 

Step 2: Set i: = i + 1 and reset i: - 1 if i > m. For i e Cj and | C } \ > l n (l n > /) examine 
whether there are clusters C p with p±j such that shifting the observation i from C } to 
C p reduces the objective function. If so, then choose C r such that the reduction 
becomes maximal and redefine Cj: = Cj — {i}, G r =C r u{i}. Otherwise return to 
step 2. 

Step 3: Repeat step 2 as long as you get any reduction, i.e. as long as i has been 
increased m times without any change. 

This stepwise optimal method works sequentially on the observations. Its result 
depends on the initial partition, on the starting observation i a and on the choice of l„. 
Normally, in order to get a suitable approximation for an optimal solution, it is 
sufficient to try several possibilities and to select the best final partition. For instance 
you can use the standard initial partition C ; = {i : i e M, i ==/ (mod n)} 5 several values 
for i 0 and /„. 

3. Implementation 

A very inefficient but numerically stable implementation of this algorithm was given 
in [5]. Using suitable up- and downdating processes for the regressions, see e.g. [1], 
'[41 results in a far more efficient algorithm. Here we have decided to use slightly 
modified FORTRAN versions of the ALGOL procedures include and regress from 
[3] working with Givens' rotations. As the downdating process may become 
numerically instable for a resulting small number of observations, see [3], it is 
recommended to use say l n &2 l 9 if m/n is large enough. As a precaution all 
internal operations are done in double precision. In order to detect numerical 
instabilities it is recommended to restart the subroutine with the found final 
partition and to compare the results. 

4. FORTRAN Subroutine 

The formal parameters of the following subroutine CWDLRS are precisely 
explained within the comment cards at the beginning of the listing. As against to 
CWDLR from [5] CWDLRS is self-contained, internally uses double precision and 
parameter communication with the central auxiliary subroutine INEXCL is done 
via COMMON statements that is faster than passing parameters via argument lists. 
Thus the two subroutines cannot directly be compared. 

SUBROUTINE CWDLKS <* f L .10 IH , N ♦ A , Y * P , KP • I A f LN »(J • X i £ t ES) 

C 

C THE PARAMETERS ARE DEFINED AS FOLLOWS • 
C 

C M NUMBER OF OBSEWVaTIO^S 

C 

C L NUMBER OF INDEPENDENT VARIABLES 

C 

C LDIH .FlkST DIMENSION OF A AND X (BELOW) IN CALLING PROGRAM 

C 

C N NUMBER OF DESIRED CLUSTERS (1 <= N <= M/LN) 

C 

C A (LDIMfN) THE ARRAY A HAS TO CONTAIN THE GIVEN <M,L)-MATRIX OF 
C OBSERVATIONS FOR Th£ INDEPENDENT VARIABLES 
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c 

C Y<M) THt A»PAY Y HAS TO CONTAIN THE GIVEN "-VECTOR OF 

C VcLUES FOR THE DEPENDENT VARIABLE 

C 

C P(M) FOR KP.NE.O THIS INTEGER M-VECTOR INITIALLY HAS TO 

C CONTAIN A FEASIBLE PARTITION OF LENGTH N VIA P(I}=J 

C (!=!♦.. .M» Jslt«,.tN). 

C FOR KP.EQ.O THE STANDARD INITIAL PARTITION IS GENERATED* 

C ON OUTPUT P WILL CONTAIN THE FINAL PARTITION OBT AINEO 

C BY THE EXCHANGE METHOD 
C 

C KP SEE P AHOVE 
C 

C tt THE EXCHANGE METHOD IS STARTED WITH OBSERVATION 

C N'JMBER IA ♦ 1 (MOO M). NORMALLY ONE SETS IA=0. 

C CHANGING I A AND USING THE SAME PART ITION P GIVES ANOTHER 

C ( BETTER ♦ EUUAL* OR WORSE) VALUE FOR THE OBJECTIVE 

C FUNCTION 

C 

C LN The MINIMUM NUMBER OF OBSERVATIONS DESIRED IN EACH 

C CLUSTER. WE MUST HAVE AT LEAST LN >= L. FOR LN <= 0 

C THIS VALUE IS AUTOMATICALLY GENERATED. IF THERE ARE 

C ENOUGH OBSERVATIONS IN PELATION^ TO THE NUMBER OF 

C CLUSTERS IT IS RECOMMENDED (ALSO FOR IMPROVING 

C NUMERICAL STABILITY) TU USE LN » L 
C 

C W(N) THE J-TH COMPONENT OF THIS INTEGER VECTOR WILL CONTAIN 

C THE NUMBER OF OBSERVATIONS I.N THE J-TH CLUSTER 

C 

C X (LDlMf N) w ILL CONTAIN THE (N,L) -MATRIX OF SOLUTION PARAMETERS, 

C- I.E. X(K»J) (K=l....tL> ARE THE REGRESSION COEFFICIENTS 

C FOR THE J-TH CLUSTER OF OBSERVATIONS (Jslr....N> 

C 

C E(N) THE J-TH COMPONENT WILL CONTAIN THE ERROR SUM OF 

C SQUARES. FOR THE J-TH CLUSTER 

C 

C £S WILL CONTAIN THE SUM OF THE E(J) 

c 
c 

C OTHER ARRAYS ARE FOR WORKING SPACE. COMMUNICATION WITH THE 

C SUBROUTINE: INEXCL IS DONE! VlA ; THE LABELED COMMON /INEX/. 

C FOR L > 10 AND/OP N > 20 DIMENSIONS HAVE TO BE ADAPTED IN 

C CWDLRS AND INEXCL •■ THE RIGHT NUMBERS ARE INDICATED IN THE 

C FOLLOWING C-CARDS 
C 

INTEGER P 1 0 « U » V ♦ W » P I tU 1 

OIMENSION A (LOlM»M) » Y (M) • X (LDIMtN) t£ <N) «P (if) tO(N) 

C: 

C LDIM <= 10 t L <= LDIM. N <= 2u 
C 

C NR= ( <L-l>*L>/2 

C 

C DIMENSION EO(N) ».XI (L).tD(LfN) »T(L»N) »R(NRtN) tOC(Ll tTC(L) * 

C • RC-(NH) tOACL) tTA<L) .RA(NR) tDfl(L) tTB(L) ?RB (NR) 

C 

REAL*8 ED (20) tXI (10) »D(10»20> »T (10t 20) »R (45* 10) tDA( 10) »f A (10) t 

* RA(45) •bB(lO) •TBI 10) tRB(45) t DC (10 ) ♦ TC ( 10 ) »RC ( 45) V 

* YI*Wl.SS»DSUMYSAtSB,SF,FF*EAfZERO.ONE*BIG 

C 

COMMON /INEX/ XI,D»T,R,6CfTC,RC»YI f WI,SS»ZER0.0NE»LLrL2»K»LU,NRV 

C 

C BIG LARGEST NUMBER ON YOUR COMPUTER 
C 

C 



C 



RIG=1.D50 

7ERO=0.D0 
ONE=l.D0 

LL=L 

L1=LL*1 

L2=LL*LL 

IF (LN.LE.O) LN=LL 
N»=( (LL-l)*LL)/2 
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C GENERATION OF INITIAL PARTITION IF DESIRED 

C 

IF(KP.NE.O) 60T0 I 
K = 0 

no 1 1=1, k 

K=K*1 

IF(K.GT.N> K=K-N 
P(I>=K 

1 CONTINUE 

C 

C INITIALIZATION TO ZERO 

C 

2 DSLM=ZERO 
DO 5 K=1.N 

g(K)=0 
EO(Kj=ZEFd 
00 3 U=I.L 

n(ijtK>=7.E«o 

T<U.<>=ZE«0 

3 CONTINUE 

•)0 4 V^lfhW 

R(V»K)=ZEtfO 
* CONTINUE 
b CONTINUE 

C 

C UPDATE FOrf INITIAL PARTITION 

C 

DO 9 I = ltf4 
K=P(I) 

IF(KP.NE.O.ANO.K.LT.O.OR.K.GT.N) RETURN 
Q(K)eQ(K)*l 

Wl=ONE 
Yl=Y<n 
DO ft U=1,L 

XI (U)sA(Utl) 

6 CONTINUE 
SS=EOCK> 
CALL INEXCL 
ED(K)sSS 

DO 7 U=ltLU 

D(U»K)=DC(U> 
T(U,K)=TC(U) 

7 CONTINUE 

DO 8 V=3 *NRV 

R<V»K)=RC<V) 

8 CONTINUE 

9 CONTINUE 

DO 10 K=ltN 

IF(0(K) .LT.LNI RETURN 
DSUMrDSUH*ED(K) 

10 CONTINUE 
IF(N.EQ.l) GOTO 22 

C 

C START OF THE EXCHANGE METhOO 

C 

I5=IA 
IT=0 

11 IS=IS*1 
IFUS.GT.M) IS=IS-m 
IF(IT.EQ.M) GOTO 22 
J=POS) 

C 

C IF THE NUMBER OF ELEMENTS OF THIS CLUSTER IS TOO SMALl 

C THEN DO NOT PEMOvE THE OBSEKVATION 

C 

IF(0(U) .L8.LN) GOTO 11 

SF=8IG 

DO 1H K=1,N 

SS=ED(K) 

Yl=Y(ISl 

00 12 11=1 fL 

XI <U)=A(U»IS) 
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12 CO^TINUF 

IFCK.EO.J) *I = - ONE 
CALL INEXCL 
IF(K.NE.J) GOTO IS 
SA=SS 

00 13 U=1.L 

OA (U) =DC (U) 
TA (U)=TC(U> 

13 CONTINUE 

DO 14 Vslv-MW 

Ra (V)=RC(V> 

14 CONTINUE 
GOTO 1H 

15 SB=i55 
FF=SB-EOtM 

IF (FF.GT.SF) GOTO ,1* 
SF=FF 

00 16 U-1»L 

OB(U)=DC(U) 
TR(U)sTC(U) 

16 CONTINUE 

00 17 V=ltNR 

Rfc(Vj=RC<V> 

17 CONTINUE 

18 CONTINUE 
EA=ED<J)-SA 
IF(5F.LT.EA) GOTO 19 

C 

C DO NOT EXCHANGE 

C 

ITslT*l 
GOTO 11 

19 IT=0 

C 

C DO EXCHANGE 

C 

P(TS)=W 

0(W>=Q(W>*1 
ED (J) =SA 
ED<W>=ED(W) >SF 
DSUM=DSUm*SF-FA 
00 2a U=1»L 

D(U*J)sDA(U) 

D(U,W)sD8(U) 

T(UrJ)sTA(U) 

T(U,to)=TB(U) 

20 CONTINUE 

00 21 V=ltNR 

RCV*J)=RA(V> 

21 CONTINUE 
GOTO 11 

C 

C CALCULATION- OF REGRESSION COEFFICIENTS FOR FINAL PARTITION 

C 

22 00 2b K s l f N 
E(K)=ED(K) 
DO 25 H=l,L 

U=Ll-tf 
TA(U)s:T(UfK) 
IF(U,EQ.L> iSOTO ?4 
NR=( <U-1>*(L2-U) >/2*l 
U1=U*1 

00 23 V=U1,L 

TA (U)=TA(U)-R(NHtK)«TA ( V) 
NR=NH+ 1 

23 CONTINUE 

24 X(U«K)=TA(U) 

25 CONTINUE 

26 CONTINUE 
ES=OSUM 
RETURN 
END 
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SUBROUTINE INEXCL 
INTEGER UtUl.V 

REAL*6 XI 110) tO (10*20) ♦ T ( 10 120 > »R (*5» 10 > tDCUO) tTCUO) ♦ 

* RC<45> ♦YItWifSSfZEROVONEf 

* XU,OU»WIXU»DPf HU»CB»SB»XVfRNf TN 

COMMON /INEX/ XI ti>tTtK v DCt TCtRCf Ylttal ,SStZERO»ONE»LLiLZf KtLUtNRV 
DO 3 U*1»LL 

IF(WI.EO.ZERO) GOTO 4 

XU=XI(U) 

IFUU.EG.ZERO) GOTO 3 

DU=DCU#K) 

WIXU=WI»XU 

DP=DU>VIXU*XU 

HU=ONE/nP 

CB=DU»HU 

SB=WIXU*HU 

DC<U)sb* 

IF<U.EQ.LL) GOTO 2 
Ul^U^l 

NR=( <U-1)»(L2-U> >/2*l 
00 i Y*U1«LL 

XVaXl (V) 

RNsR ( MR a K ) 

XI (V) =XV-XU»RN 

RC (NfO=CB»RN*SB*XV 

NRV=NR 

NR=NR* 1 

1 CONTINUE 

2 XV=YI 
TN=T(UtK) 
YI=XV-XU*TN 
TC(U>=CB*TN*SB*XV 
LU=U 

3 CONTINUE 
SS=5S*nI*YI#YI 

4 RETURN 
END 



5. Computing Time 

As the computing time heavily depends oft the initial partition, on i a and on l H , it is 
not possible to give a precise estimate. Normally about six passes through the 
observations, are enough. For an example with m=96, Z = 5, /„ — Z, n = 2,3,4, i fl = 0, 
KP=0 (standard initial partition) CWDLRS has needed about 30 seconds on a 
TR 440 computer (about half as fast as an IBM 370/158). In this case CWDLRS was 
about 20 times faster than the corrected version of CWDLR, see [5], [7]. The gain 
will be larger for higher values of /, n, and m. 
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